Dispersity-Driven Melting Transition in Two Dimensional Solids 
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We perform extensive simulations of 10* Lennard- Jones particles to study the effect of particle 
size dispersity on the thermodynamic stability of two-dimensional solids. We find a novel phase 
diagram in the dispersity-density parameter space. We observe that for large values of the density 
there is a threshold value of the size dispersity above which the solid melts to a liquid along a 
line of first order phase transitions. For smaller values of density, our results are consistent with 
the presence of an intermediate hexatic phase. Further, these findings support the possibility of a 
multicritical point in the dispersity-density parameter space. 
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Recently there has been considerable interest in what 
happens to the liquid-sohd transition in a system if the 
constituent particles are not all identical but have dif- 
ferent sizes. The question was first raised in the con- 
text of colloidal solutions [Q, and subsequently addressed 
for other systems These studies mainly focused 

on the effect of size dispersity A on the P — p equa- 
tion of state, where P and p denote pressure and density. 
On increasing A from zero, the density discontinuity at 
the transition decreases, eventually vanishing at a critical 
value A = Ac above which there is no liquid-solid density 
discontinuity. This remarkable phenomenon — similar to 
the effect of temperature T on the conventional liquid-gas 
phase transition|^ — occurs in both two and three dimen- 
sions, and for various forms of interaction potentials and 
size distributions [Q. 

These seminal studies leave some questions unan- 
swered. First, what are the structures of the phases? 
Second, can one pass continuously from solid to liq- 
uid "around the critical point" at Ac, just as one can 
pass continuously from liquid to gas "around the critical 
point" at Tc? A "yes" answer would not be consistent 
with the common picture of melting as a first order phase 
transition (which cannot have a critical point because of 
symmetry mismatch of the two phases Q). A "no" an- 
swer would lead to a natural third question: In the A — p 
parameter space, what is the location and nature of the 
phase boundary between crystalline and liquid phases? 
The third question has not gone unnoticed — indeed, Ref. 

simulates a binary mixture of 108 "soft" disks, and 
shows that upon increasing A the crystal undergoes a 
transition to an amorphous solid at a threshold disper- 
sity At/i, suggesting that the transition is of first order. 

Here we address all three questions by simulating a 
relatively large system comprised oi N = 10** Lennard- 
Jones particles of two different radii in a square box of 
edge Lq with periodic boundary condition. With each 
particle i, we associate a size parameter a^, and define 
the distance scale for the interaction between particles i 
and j to be aij = ai + aj. We assign to half the par- 
ticles the value Ui = cto(1 + A), and to the other half 



the value ai — aQ{l — A). If particles i and j are at a 
distance smaller than a cutoff distance Tc, they in- 
teract via a "shifted-force Lennard- Jones" potential |q| 
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(cTy/rij)^^ - (<7y/r-ijTJ + finj). Here /(r^) 
is a linear function whose coefficients are chosen such 
that and its gradient, the force, continuously van- 



ish at 



Since takes its minimum value at 
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we consider this equilibrium 



distance to be the sum of the radii of the two particles 
i and j, Rij = Ri + Rj, so the radius of particle i is 
Ri = 2^1^ X Oi and the average radius is 2^1'° x (Tq- 

We perform molecular dynamics (MD) simulations us- 
ing the velocity Verlet integrator method We record 
the results in reduced units in which ctq is 0.5, and the 
Lennard-Jones energy scale e, the particle mass, and 
Boltzmann constant are all unity. In these units, we 
choose Tc = 2.5 and the length of each MD time step 
8t — O.Ol. The system is first thermalized at T = 1, 
using the Berendsen rescaling method [||, for a period 
of length t; typically r = (5 x 10'')(5<. Then we run 
the system for an additional period r as a constant NVE 
system (micro-canonical ensemble) . We continuously cal- 
culate P, T and energy E, and we consider the system to 
be in equilibrium only when the fluctuations of all three 
quantities are less than 1% of their average values. The 
thermalization time r is chosen to be more than the time 
it takes for the system to equilibrate. 

We define the size dispersity to be the ratio of the size 
distribution variance to its average which equals A 
in our model, and we define p = )/^0' the ra- 

tio of the total area assigned to the disks to the system 
area. For each value of A, we start by placing the lO'' 
particles randomly on the sites of a square lattice of edge 
io ~ 150; higher density states are obtained by gradu- 
ally compressing the system by reducing Lq- Typically 
the starting density is p = 0.7, and we increase p to 1.05 
through approximately 10 intermediate densities, equili- 
brating the system at each[^. 

We present our results for the state points with T = 1, 
p = 0.90 - 1.05 and A = - 0.12. At these densities, 
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the A = system is a 2D-solid with a triangular or- 
der, but at large A the system becomes disordered and 
a liquid. By probing the translational and orientational 
order, we determine the phase of each state point and we 
locate the transition between the two phases. To study 
translational order, we calculate the total pair correla- 
tion function g{r), as well as the partial functions gii{r), 
922{r) and gi2{r) Here g{r) is the probabihty dis- 
tribution of finding two particles at a distance r, and 
gij{r) is the same for an pair (i = 1 stands for 

small and i = 2 for large particles). We find that all 
three gij{r) display behavior similar to indicating 
that the system maintains its substitutionally-disordered 
configuration and does not tend toward de-mixing. 

In Fig. |l|(a) we show the effect of tuning A on transla- 
tional order. We observe that the monodisperse (A = 0) 
system shows the quasi-long-range translational order ex- 
pected for a 2D-solid| 10|, characterized by a power-law 
decay of the envelope of g{r) and the persistence of the 
solid structure periodicity up to very large distances. For 
A < Ath{p), where Ath{p) is the threshold value at fixed 
p, the solid maintains this quasi-long-range order, al- 
though the decay exponent appears to increase somewhat 
with A. For A > Ath{p), we observe a qualitative change 
in the structure: the quasi-long-range translational order 
disappears, and is replaced by an exponential decay of 
the envelope of g{r), which at very long distances shows 
the uniform distribution of a structureless liquid. We ob- 
serve this behavior for all densities between p = 0.96 and 
p = 1.05, and find that 0.09 < Athip) < 0.10 for aU p. 

Next we study the local bond orientational order by 
calculating for each particle j the sixfold orientational 
order parameter Q 
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The sum runs over all z nearest neighbors k oi j, and 6jk 
is the angle of the bond joining particles j and k with 
respect to a fixed axis. We identify the nearest neigh- 
bors as the particles that are closer than the location of 
the first minimum of g{r). The modulus of will be 
unity if the neighbors form a perfect hexagon around j, 
which occurs for all particles in a triangular lattice, the 
close-packed configuration of a 2D-solid. For a distorted 
hexagon or a different polygon, < 1 — e.g. for a 

liquid, the distribution of \{'4')j \ centers around 0.5 Q. 

We define the continuous order parameter field ■!/'(r)as 
the value of if the position of particle j is Vj = r, and 
we calculate the orientational correlation functionlllSl 



ge{\r - ro|) = {■ip{r)^p{ro)), 
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where (...) denotes an average over r,ro and time. 
Fig. |l|(b) shows that if A is small and p is large, the sys- 
tem displays the long-range orientational order of a solid 
in that luiir—,00 gei''') 7^ 0- Noteworthy is that for each 
value of p, orientational order disappears upon a small in- 
crease in dispersity near Ath{p)- For A > Ath{p), g&ir) 



appears to decay exponentially, which identifies the sys- 
tem as a liquid. Similar plots for other large values of 
p suggest that near p — 1.0, there is a first order phase 
transition from solid to liquid, driven by an increase in 
A. This observation is in agreement with the results of 
Ref. §. Fig.! shows that the small dispersity system has 
the ordered structure of a solid while the large dispersity 
system has the disordered structure of a liquid — providing 
an answer to the first of the three questions. 

Since identifying phases relies on the behavior of the 
system in the thermodynamic limit, we apply finite size 
scaling to the moments of the orientational order param- 
eter if], which is the average value of We use stan- 
dard techniques originally developed for the Ising model 
1 14 1, and recently applied to the 2D melting transition 
|15,|l^. In order to calculate the moments of -0 at a 
scale b = Lq/M, we divide the system into blocks of 
edge b and we define ipb for each block as the absolute 
value of the average of ip{v) over the block. Then we find 
the moments of 4'b by averaging over all blocks and all 
configurations of the system after equilibration jl^]. 

To explore the precise location of the phase transition, 
we calculate the cumulant [|l^ 
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For a completely ordered (solid) system, Ub — 2/3 in the 
thermodynamic limit, while for a disordered (liquid) sys- 
tem Ub — > 0. For an infinite system, Ub jumps between 
these two limiting values at the phase transition point 
and for finite systems, this jump becomes rounded. Still, 
one can determine the location of a phase transition by 
finding the point at which the Ub curves for different sys- 
tem sizes intersect ||l5[. In Fig. ||, we plot Ub versus A, 
for different scales 6 at a fixed p. We find the transition 
from the value 2/3 to lower values upon passing through 
the phase transition. Moreover, we estimate Ath{p) from 
the crossing point of the curves. In the phase diagram 
of Fig. ^, the line Lth is the locus of all such threshold 
points separating solid and liquid phases and shows that 
for p > 0.96, Ath{p) ~ 0.097 is essentially independent 
of p. 

Next we study the finite size scaling of (ipb)- Because 
of the qualitative difference in the form of g^ (r) between 
solid and liquid, the behavior of (^ f ) a s a function of 
b changes drastically upon melting |16|. In the liquid 



for 6 ^ where ^ is the correlation length, {ip'^) de- 
cays as b^^, while for the solid, {tp'^) remains constant. 
Fig. ^(a) shows that for p = 1.0, the behavior of the 
system changes abruptly from solid (given by the line 
with zero slope) to liquid at Ath, which is consistent 
with our previous studies of g{r), geir) and Ub- The 
A = 0.10 liquid curve shows long-range correlation for 
b < which crosses over to short-range correlation (slope 
—2) for 6 » ^ (^ « O.6L0 for this curve). The correspond- 
ing plots for larger A show that ^ shrinks upon increasing 
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A. For p = 0.9 (Fig. |(b)), we observe both solid behav- 
ior for A < 0.06 and Hquid behavior for A > 0.06. For 
A = 0.06, Fig. ^(b) shows an algebraic decay for the cor- 
relation function, with exponent —1/4. This "intermedi- 
ate" behavior is reminiscent of the hexatic phase for 
which the orientational correlation decays algebraically 
while the system does not possess quasi-long-range trans- 
lational order. In Fig. ^ we have specified as the / phase 
the locus of the points showing this intermediate behav- 
ior. 

In summary, we have studied a melting transition 
driven not by T but by A. We have simulated rela- 
tively large systems and applied finite size scaling (Figs ||, 
^ arriving at a phase diagram for this dispersity-driven 
melting (Fig. ^). Melting takes place from a 2D or- 
dered phase to a disordered liquid phase, similar to 
the conventional temperature-driven melting processes. 
Moreover, at large values of p, melting is a first order 
phase transition at a threshold dispersity value Ath{p) = 
0.097±0.005. Our study of the mean square displacement 
of the particles shows that this melting is accompanied 
by a transition from a frozen solid to a diffusive liquid, 
distinguishing it from the glass transition observed in . 
The threshold line Lth extends almost horizontally down 
to the point C with coordinates Ac ~ 0.097, pc ~ 0.96. 
Below pc, finite size scaling of the orientational order pa- 
rameter suggests the existence of an intermediate "hex- 
atic" phase between the solid and liquid phases. Thus 
we hypothesize that point C is a multicritical point, 
where two lines of continuous transitions (separating 
liquid/ "hexatic" and "hexatic" /solid phases) meet the 
line of first order transitions Lth (separating liquid/solid 
phases) as shown in Fig. Fig. || provides an an- 

swer to the last two of the three questions: one can not 
pass continuously from solid to liquid "around the crit- 
ical point" C, because the two phases are separated by 
the line of first order phase transitions Lth- A similar 
horizontal line of order-disorder transition has been ob- 
served in the study of the effect of quenched impurities 
on the structure of 2D-solids jl^ . 
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tremely helpful criticism, NSF for financial support and 
the Boston University Center for Computational Science 
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FIG. 1. EfTect of dispersity A on tramslational and orien- 
tational order at p = 1.0. (a)Total pair correlation function 
g{r) as a function of distance r. All curves oscillate around the 

value g{'r) = 1, so wc have separated them to facilitate com- 
parison. We find that a transition from solid to liquid struc- 
ture occurs on increasing the dispersity between A = 0.09 and 
A — 0.10. (b) Normalized oriontational correlation function 
versus r. The change between A = 0.09 and A = 0.10 cor- 
responds to the transition from an orientationally long-range 
correlated solid to a short-range correlated liquid. Both sets 
of curves show that for p = 1.0 the solid-liquid transition oc- 
curs at a value of dispersity between 0.09 and 0.10 which is 
consistent with the observations based on Figs. 2 — 4. 
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FIG. 2. Cumulant Ut of the bond oriontational order pa- 
rameter ^ as a function of dispersity A for p = 1.0. Different 
curves correspond to different scales 6, where b = Lq/M is 
the block size, so smaller M corresponds to larger scale. The 
dotted lines connecting the data points are guides to the eye. 
We identify the threshold value of A to be the point where 
all the curves for different scales intersect. 
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FIG. 3. Phase diagram in the A — p (dispersity-donsity) 
parameter space. Squares represent solid points and circles 
liquid points. The threshold line Lth connects crosses, which 
are the first order phase transition points derived from the 
cumulant analysis. The triangles are the points of the inter- 
mediate (/) phase, showing a hexatic behavior. The large 
diamond marks the multicritical point C. 
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FIG. 4. Log-log plots of (i/Jb) versus 6 for different disper- 

sities A. Both axes are normalized by their values at system 
size Lq, which causes the curves to meet at the origin and fa- 
cilitates comparing their asymptotic slopes. Dashed lines con- 
necting the data points are guides to the eye. Solid straight 
lines arc reference lines with slopes —2 and —1/4. (a) For 
p = 1.0 there is an abrupt change from asymptotic slope of 
(data lying on the abscissa) to —2, which corresponds to a 
solid-liquid transition on increasing the dispersity above At/j. 
(b) For p = 0.9, the A = 0.06 curve falls between the solid 
and liquid regimes, and shows a slope of —1/4 characteristic 
of the hexatic phase. 
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